//! 多项式的计算

const std = @import("std");
const lbmath = @import("./lbmath.zig");
const Allocator = std.mem.Allocator;
const dprint = std.debug.print;
test "for loop" {
    for (0..3) |i| {
        dprint(" for  loop:{}\n", .{i});
    }
}
pub const Polynomial_f64 = Polynomial(f64);

pub fn Polynomial(comptime T: type) type {
    return struct {
        allocator: Allocator,
        const Self = @This();
        // 结构体内部的test函数默认不执行
        test "polynomial test" {
            dprint("polynomial test************\n", .{});
        }
        // test addpolynomials {
        //     var arena = std.heap.ArenaAllocator.init(std.heap.page_allocator);
        //     defer arena.deinit();
        //     const allocator = arena.allocator();
        //     const polynomial_f64 = Polynomial(f64).init(allocator);
        // }
        pub fn init(allocator: Allocator) Self {
            return Self{
                .allocator = allocator,
            };
        }
        // 多项式相加
        // void addpolynomials(Polynomial* p1, Polynomial* p2,Polynomial* result)
        pub fn addpolynomials(self: Self, p1: []T, p2: []T, result: ?[]T) ?[]T {
            const maxDegree = if (p1.len > p2.len) p1.len else p2.len;
            const re = if (result) |re1|
                if (re1.len < maxDegree) return null else re1
            else
                self.allocator.alloc(T, maxDegree) catch return null;
            for (0..maxDegree) |i| {
                re[i] = (if (i < p1.len) p1[i] else 0) + (if (i < p2.len) p2[i] else 0);
            }
            return re;
        }

        // 多项式相乘
        pub fn multiplypolynomials(self: Self, p1: []T, p2: []T, result: ?[]T) ?[]T {
            const resultDegree = p1.len + p2.len;
            const re = if (result) |re1|
                if (re1.len < resultDegree) return null else re1
            else
                self.allocator.alloc(T, resultDegree) catch return null;

            @memset(re, 0);
            for (0..p1.len) |i| {
                for (0..p2.len) |j| {
                    re[i + j] += p1[i] * p2[j];
                }
            }
            return re;
        }

        fn polynomial_plus_pow_xan(p: []T, a: T, fac: T, degree: u32) void {
            // 从degree中选
            // 0 .. degree
            var temp: T = 1;
            var cnm: u32 = 1;
            var i: i32 = degree;
            while (i >= 0) : (i -= 1) {
                p[i] += (fac * temp * cnm);
                temp *= a;
                const j: u32 = degree - i + 1;
                cnm *= ((degree - j + 1) / j);
            }
        }

        // x根正方向移动是 x-offset
        // 负方向移动是 x+offset
        //
        pub fn polynomialshift_root(self: Self, p: []T, offset: T, result: ?[]T) ?[]T {
            const re: []T = if (result) |re1|
                if (re1.len < p.len) return null else re1
            else
                self.allocator.alloc(T, p.len) catch return null;
            @memset(re, 0);
            for (0..re.len) |i| {
                polynomial_plus_pow_xan(re, -offset, p[i], i);
                // p->coefficients[i]*=fac;
            }

            return re;
        }
        pub fn scala_mult(self: Self, p: []T, fac: T, result: ?[]T) ?[]T {
            const re: []T = if (result) |re1|
                if (re1.len < p.len) return null else re1
            else
                self.allocator.alloc(T, p.len) catch return null;
            for (0..p.len) |i| {
                re[i] = p[i] * fac;
            }
            return re;
        }
        // 求解一次方程: a*x + b = 0
        pub fn solve_linear(self: Self, coeffs: []T, roots: ?[][2]T) ?[][2]T {
            if (coeffs.len < 2) return null;
            const a = coeffs[0];
            const b = coeffs[1];

            if (b == 0) return null; // 不是一次方程

            const re: [][2]T = if (roots) |re1|
                if (re1.len < 1) return null else re1
            else
                self.allocator.alloc([2]T, 1) catch return null;

            re[0][0] = -a / b;
            re[0][1] = 0;
            return re;
        }
        pub fn solve_quadratic(self: Self, coeffs: []T, roots: ?[][2]T) ?[][2]T {
            if (coeffs.len < 3) return null;
            const c = coeffs[0];
            const b = coeffs[1];
            const a = coeffs[2];

            if (a == 0) return solve_linear(self, coeffs, roots);
            const disc: T = b * b - 4 * a * c;

            const re: [][2]T = if (roots) |re1|
                if (re1.len < 2) return null else re1
            else
                self.allocator.alloc([2]T, 2) catch return null;

            if (disc >= 0) {
                const sqrt_disc = std.math.sqrt(disc);
                re[0][0] = (-b - sqrt_disc) / (2 * a);
                re[0][1] = 0;
                re[1][0] = (-b + sqrt_disc) / (2 * a);
                re[1][1] = 0;
            } else {
                const real_part = -b / (2 * a);
                const imag_part = std.math.sqrt(-disc) / (2 * a);
                re[0][0] = real_part;
                re[0][1] = imag_part;
                re[1][0] = real_part;
                re[1][1] = -imag_part;
            }

            return re;
        }
        pub fn solve_cubic(self: Self, coeffs: []T, roots: ?[][2]T) ?[][2]T {
            if (coeffs.len < 4) return null;
            const d = coeffs[0];
            const c = coeffs[1];
            const b = coeffs[2];
            const a = coeffs[3];

            if (std.math.approxEqAbs(T, a, 0, null)) return self.solve_quadratic(coeffs[0..3], roots);

            const re: [][2]T = if (roots) |re1|
                if (re1.len < 3) return null else re1
            else
                self.allocator.alloc([2]T, 3) catch return null;

            // 归一化
            const p = (3 * a * c - b * b) / (3 * a * a);
            const q = (2 * b * b * b - 9 * a * b * c + 27 * a * a * d) / (27 * a * a * a);

            const discriminant = (q * q) / 4 + (p * p * p) / 27;

            // 使用容差判断判别式
            const eps = lbmath.defaultTolerance(T);
            if (discriminant > eps) {
                // 一个实根，两个共轭复根
                const sqrt_disc = std.math.sqrt(discriminant);
                const u_cubed = -q / 2 + sqrt_disc;
                const v_cubed = -q / 2 - sqrt_disc;

                const u = std.math.cbrt(u_cubed);
                const v = std.math.cbrt(v_cubed);

                const x1 = u + v - b / (3 * a);
                re[0][0] = x1;
                re[0][1] = 0;

                // 复根
                const real_part = -(u + v) / 2 - b / (3 * a);
                const imag_part = (std.math.sqrt(3) / 2) * (u - v);
                re[1][0] = real_part;
                re[1][1] = imag_part;
                re[2][0] = real_part;
                re[2][1] = -imag_part;
            } else if (discriminant < -eps) {
                // 三个实根
                // 确保参数非负（处理浮点误差）
                const safe_p = -@abs(p);
                const r = std.math.sqrt(@abs(safe_p * safe_p * safe_p / 27));

                // 确保 acos 的参数在 [-1, 1] 范围内
                const phi_arg = -q / (2 * r);
                const phi = std.math.acos(phi_arg);

                const s = 2 * std.math.sqrt(@abs(safe_p / 3));

                re[0][0] = s * std.math.cos(phi / 3) - b / (3 * a);
                re[0][1] = 0;
                re[1][0] = s * std.math.cos((phi + 2 * std.math.pi) / 3) - b / (3 * a);
                re[1][1] = 0;
                re[2][0] = s * std.math.cos((phi + 4 * std.math.pi) / 3) - b / (3 * a);
                re[2][1] = 0;
            } else {
                // 判别式 ≈ 0，有重根
                if (std.math.approxEqAbs(T, q, 0, null)) {
                    // q ≈ 0，p ≈ 0，三重实根
                    const root = -b / (3 * a);
                    re[0][0] = root;
                    re[0][1] = 0;
                    re[1][0] = root;
                    re[1][1] = 0;
                    re[2][0] = root;
                    re[2][1] = 0;
                } else {
                    // 一个二重根和一个单根
                    const u = std.math.cbrt(-q / 2);
                    const x1 = 2 * u - b / (3 * a);
                    const x2 = -u - b / (3 * a);

                    re[0][0] = x1;
                    re[0][1] = 0;
                    re[1][0] = x2;
                    re[1][1] = 0;
                    re[2][0] = x2;
                    re[2][1] = 0;
                }
            }
            return re;
        }
        pub fn solve_quartic(self: Self, coeffs: []T, roots: ?[][2]T) ?[][2]T {
            const eps = lbmath.defaultTolerance(T); // 统一误差阈值

            if (coeffs.len < 5) return null;
            const e = coeffs[0];
            const d = coeffs[1];
            const c = coeffs[2];
            const b = coeffs[3];
            const a = coeffs[4];

            if (std.math.approxEqAbs(T, a, 0, null))
                return self.solve_cubic(coeffs[0..4], roots);

            const re: [][2]T = if (roots) |re1|
                if (re1.len < 4) return null else re1
            else
                self.allocator.alloc([2]T, 4) catch return null;

            // 归一化系数 a*x^4 + b*x^3 + c*x^2 + d*x + e = 0
            const a_norm = a;
            const b_norm = b / a_norm;
            const c_norm = c / a_norm;
            const d_norm = d / a_norm;
            const e_norm = e / a_norm;

            // 通过替换 x = y - b/4 转换为 depressed quartic: y^4 + p*y^2 + q*y + r = 0
            const b_over_4 = b_norm / 4;
            const b_sq = b_norm * b_norm;
            const b_cube = b_sq * b_norm;
            const b_quad = b_cube * b_norm;

            const p = c_norm - 6 * b_sq / 16;
            const q = d_norm - b_norm * c_norm / 2 + b_cube / 8;
            const r = e_norm - b_norm * d_norm / 4 + b_sq * c_norm / 16 - 3 * b_quad / 256;

            // 使用 Ferrari 方法
            if (std.math.approxEqAbs(T, q, 0, null)) {
                // 双二次方程: y^4 + p*y^2 + r = 0
                var quad_coeffs: [3]T = .{ r, p, 1 }; // 注意系数顺序: [r, p, 1] 对应 z^2 + p*z + r
                var quad_roots1: [2][2]T = .{[2]T{ 0, 0 }} ** 2;

                const quad_roots = self.solve_quadratic(&quad_coeffs, &quad_roots1) orelse {
                    @memset(re, .{ 0, 0 });
                    return re;
                };

                var root_count: usize = 0;
                for (quad_roots) |z_root| {
                    const z_real = z_root[0];
                    const z_imag = z_root[1];

                    if (std.math.approxEqAbs(T, z_imag, 0, null)) {
                        // z 是实数
                        if (z_real >= -eps) { // 允许微小负值
                            const abs_z = if (z_real < 0) -z_real else z_real;
                            const sqrt_z = std.math.sqrt(abs_z);

                            // 正根
                            if (z_real >= 0) {
                                const real_sqrt_z = std.math.sqrt(z_real);
                                if (root_count < 4) {
                                    re[root_count] = .{ real_sqrt_z - b_over_4, 0 };
                                    root_count += 1;
                                }
                                // ✅ 修复1：始终添加第二个根（即使 z=0）
                                if (root_count < 4) {
                                    re[root_count] = .{ -real_sqrt_z - b_over_4, 0 };
                                    root_count += 1;
                                }
                            }
                            // z < 0 → 复数根
                            else if (root_count < 4) {
                                re[root_count] = .{ -b_over_4, sqrt_z }; // +i√|z|
                                root_count += 1;
                                if (root_count < 4) {
                                    re[root_count] = .{ -b_over_4, -sqrt_z }; // -i√|z|
                                    root_count += 1;
                                }
                            }
                        }
                    } else {
                        // z 是复数：使用稳定复数开方算法
                        const rr = std.math.hypot(T, z_real, z_imag); // |z|
                        if (rr < eps) {
                            // ✅ 修复2：添加两个相同的根
                            if (root_count < 4) {
                                re[root_count] = .{ -b_over_4, 0 };
                                root_count += 1;
                            }
                            if (root_count < 4) {
                                re[root_count] = .{ -b_over_4, 0 };
                                root_count += 1;
                            }
                            continue;
                        }

                        // 复数开方公式
                        const sqrt_r = std.math.sqrt(rr);
                        const theta = std.math.atan2(T, z_imag, z_real);
                        const half_theta = theta / 2;
                        const real_part = sqrt_r * @cos(half_theta);
                        const imag_part = sqrt_r * @sin(half_theta);

                        // 两个平方根: ±(real_part + i*imag_part)
                        if (root_count < 4) {
                            re[root_count] = .{ real_part - b_over_4, imag_part };
                            root_count += 1;
                        }
                        if (root_count < 4) {
                            re[root_count] = .{ -real_part - b_over_4, -imag_part };
                            root_count += 1;
                        }
                    }
                }

                // 填充剩余根为 (0,0)
                while (root_count < 4) : (root_count += 1) {
                    re[root_count] = .{ 0, 0 };
                }
            } else {
                // 一般情况：构造辅助三次方程: m^3 + (2*p)*m^2 + (p^2 - 4*r)*m - q^2 = 0
                var cubic_coeffs = [4]T{ 1, 2 * p, p * p - 4 * r, -q * q };
                var cubic_roots1: [3][2]T = .{[2]T{ 0, 0 }} ** 3;
                const cubic_roots = self.solve_cubic(&cubic_coeffs, &cubic_roots1) orelse {
                    @memset(re, .{ 0, 0 });
                    return re;
                };

                // 选择最佳实根 m (优先选择使 s_squared > 0 的根)
                var m: T = 0;
                var found_m = false;
                var best_s_sq: T = -std.math.inf(T);

                for (cubic_roots) |root| {
                    if (std.math.approxEqAbs(T, root[1], 0, null)) { // 实根
                        const m_candidate = root[0];
                        const s_sq_candidate = 2 * m_candidate - p;

                        if (s_sq_candidate > 0) {
                            // 选择最大的 s_sq_candidate
                            if (!found_m || s_sq_candidate > best_s_sq) {
                                m = m_candidate;
                                best_s_sq = s_sq_candidate;
                                found_m = true;
                                // 优先选择明显为正的根
                                if (s_sq_candidate > eps) break;
                            }
                        }
                    }
                }

                if (!found_m || best_s_sq < -eps) {
                    // 无有效实根，返回空解
                    @memset(re, .{ 0, 0 });
                    return re;
                }

                // 计算分解参数（处理 s ≈ 0 的情况）
                const s_squared = 2 * m - p;
                if (s_squared < eps) {
                    // ✅ 修复3：允许 s_squared > 0 但 < eps 的情况
                    if (s_squared <= -eps) {
                        // s_squared 为负且超出容差，返回空解
                        @memset(re, .{ 0, 0 });
                        return re;
                    }
                    // 否则继续，s_squared 很小但为正
                }

                const s_val = std.math.sqrt(s_squared);
                const q_scale = if (s_val > eps) q / s_val else 0;
                const t1 = (s_squared + p - q_scale) / 2;
                const t2 = (s_squared + p + q_scale) / 2;

                // 解两个二次方程
                var root_count: usize = 0;

                // 第一个二次方程: y^2 + s_val*y + t1 = 0
                const disc1 = s_val * s_val - 4 * t1;
                if (disc1 >= -eps) {
                    const sqrt_d = if (disc1 > eps)
                        std.math.sqrt(disc1)
                    else
                        0; // 处理微小负判别式

                    // 始终添加两个根（即使重根）
                    if (root_count < 4) {
                        re[root_count] = .{ (-s_val + sqrt_d) / 2 - b_over_4, 0 };
                        root_count += 1;
                    }
                    if (root_count < 4) {
                        re[root_count] = .{ (-s_val - sqrt_d) / 2 - b_over_4, 0 };
                        root_count += 1;
                    }
                } else {
                    // 复根
                    const sqrt_d = std.math.sqrt(-disc1);
                    if (root_count < 4) {
                        re[root_count] = .{ -s_val / 2 - b_over_4, sqrt_d / 2 };
                        root_count += 1;
                    }
                    if (root_count < 4) {
                        re[root_count] = .{ -s_val / 2 - b_over_4, -sqrt_d / 2 };
                        root_count += 1;
                    }
                }

                // 第二个二次方程: y^2 - s_val*y + t2 = 0
                const disc2 = s_val * s_val - 4 * t2;
                if (disc2 >= -eps) {
                    const sqrt_d = if (disc2 > eps)
                        std.math.sqrt(disc2)
                    else
                        0; // 处理微小负判别式

                    if (root_count < 4) {
                        re[root_count] = .{ (s_val + sqrt_d) / 2 - b_over_4, 0 };
                        root_count += 1;
                    }
                    if (root_count < 4) {
                        re[root_count] = .{ (s_val - sqrt_d) / 2 - b_over_4, 0 };
                        root_count += 1;
                    }
                } else {
                    // 复根
                    const sqrt_d = std.math.sqrt(-disc2);
                    if (root_count < 4) {
                        re[root_count] = .{ s_val / 2 - b_over_4, sqrt_d / 2 };
                        root_count += 1;
                    }
                    if (root_count < 4) {
                        re[root_count] = .{ s_val / 2 - b_over_4, -sqrt_d / 2 };
                        root_count += 1;
                    }
                }

                // 用0填充剩余位置
                while (root_count < 4) : (root_count += 1) {
                    re[root_count] = .{ 0, 0 };
                }
            }

            return re;
        }
        // pub fn solve_quartic(self: Self, coeffs: []T, roots: ?[][2]T) ?[][2]T {
        //     if (coeffs.len < 5) return null;
        //     const e = coeffs[0];
        //     const d = coeffs[1];
        //     const c = coeffs[2];
        //     const b = coeffs[3];
        //     const a = coeffs[4];

        //     if (a == 0) return self.solve_cubic(coeffs[0..4], roots);

        //     const re: [][2]T = if (roots) |re1|
        //         if (re1.len < 4) return null else re1
        //     else
        //         self.allocator.alloc([2]T, 4) catch return null;

        //     // 归一化系数 a*x^4 + b*x^3 + c*x^2 + d*x + e = 0
        //     const a_norm = a;
        //     const b_norm = b / a_norm;
        //     const c_norm = c / a_norm;
        //     const d_norm = d / a_norm;
        //     const e_norm = e / a_norm;

        //     // 通过替换 x = y - b/4 转换为 depressed quartic: y^4 + p*y^2 + q*y + r = 0
        //     const b_over_4 = b_norm / 4;
        //     const b_sq = b_norm * b_norm;
        //     const b_cube = b_sq * b_norm;
        //     const b_quad = b_cube * b_norm;

        //     const p = c_norm - 6 * b_sq / 16;
        //     const q = d_norm - b_norm * c_norm / 2 + b_cube / 8;
        //     const r = e_norm - b_norm * d_norm / 4 + b_sq * c_norm / 16 - 3 * b_quad / 256;

        //     // 使用 Ferrari 方法
        //     if (std.math.approxEqAbs(T, q, 0, 1e-12)) {
        //         // 双二次方程: y^4 + p*y^2 + r = 0
        //         // 双二次方程: y^4 + p*y^2 + r = 0
        //         var quad_coeffs: [3]T = [_]T{ 1, p, r };
        //         var quad_roots1: [2][2]T = [_][2]T{[2]T{ 0, 0 }} ** 2;

        //         const quad_roots = self.solve_quadratic(&quad_coeffs, &quad_roots1) orelse {
        //             @memset(@as([*][2]T, @ptrCast(re))[0..4], [2]T{ 0, 0 });
        //             return re;
        //         };

        //         var root_count: usize = 0;
        //         for (quad_roots) |z_root| {
        //             const z_real = z_root[0];
        //             const z_imag = z_root[1];

        //             if (std.math.approxEqAbs(T, z_imag, 0, 1e-12)) {
        //                 // z 是实数
        //                 if (z_real >= 0) {
        //                     const sqrt_z = std.math.sqrt(z_real);
        //                     // 正根
        //                     if (root_count < 4) {
        //                         re[root_count][0] = sqrt_z - b_over_4;
        //                         re[root_count][1] = 0;
        //                         root_count += 1;
        //                     }
        //                     // 负根（如果非零）
        //                     if (z_real > 0 and root_count < 4) {
        //                         re[root_count][0] = -sqrt_z - b_over_4;
        //                         re[root_count][1] = 0;
        //                         root_count += 1;
        //                     }
        //                 } else {
        //                     // z < 0 → y = ±i√|z|，实部为0，虚部为 ±√|z|
        //                     const sqrt_neg_z = std.math.sqrt(-z_real);
        //                     if (root_count < 4) {
        //                         re[root_count][0] = -b_over_4; // 实部
        //                         re[root_count][1] = sqrt_neg_z; // 虚部
        //                         root_count += 1;
        //                     }
        //                     if (root_count < 4) {
        //                         re[root_count][0] = -b_over_4; // 实部
        //                         re[root_count][1] = -sqrt_neg_z; // 虚部
        //                         root_count += 1;
        //                     }
        //                 }
        //             } else {
        //                 // z 是复数：解 y^2 = z（复数开方）
        //                 const rr = std.math.sqrt(z_real * z_real + z_imag * z_imag); // |z|
        //                 const real_part = std.math.sqrt((rr + z_real) / 2); // Re(√z)
        //                 const imag_part = std.math.sqrt((rr - z_real) / 2); // Im(√z)
        //                 const sign = if (z_imag >= 0) T(1) else T(-1); // 符号

        //                 // √z = ± (real_part + i * sign * imag_part)
        //                 if (root_count < 4) {
        //                     re[root_count][0] = real_part - b_over_4; // Re(y)
        //                     re[root_count][1] = sign * imag_part; // Im(y)
        //                     root_count += 1;
        //                 }
        //                 if (root_count < 4) {
        //                     re[root_count][0] = -real_part - b_over_4; // Re(y)
        //                     re[root_count][1] = -sign * imag_part; // Im(y)
        //                     root_count += 1;
        //                 }
        //             }
        //         }

        //         // 填充剩余根为 (0,0)
        //         while (root_count < 4) : (root_count += 1) {
        //             re[root_count][0] = 0;
        //             re[root_count][1] = 0;
        //         }
        //     } else {
        //         // 一般情况：使用 Ferrari 方法
        //         // 构造辅助三次方程: m^3 + (2*p)*m^2 + (p^2 - 4*r)*m - q^2 = 0

        //         var cubic_coeffs = [_]T{0} ** 4;
        //         cubic_coeffs[0] = 1;
        //         cubic_coeffs[1] = 2 * p;
        //         cubic_coeffs[2] = p * p - 4 * r;
        //         cubic_coeffs[3] = -q * q;
        //         var cubic_roots1 = [_][2]T{[2]T{ 0, 0 }} ** 3;
        //         const cubic_roots = self.solve_cubic(cubic_coeffs, &cubic_roots1) orelse {
        //             @memset(@as([*][2]T, @ptrCast(re))[0..4], [2]T{ 0, 0 });
        //             return re;
        //         };

        //         // 取第一个实根作为 m
        //         var m: T = undefined;
        //         var found_m = false;
        //         for (cubic_roots) |root| {
        //             if (std.math.approxEqAbs(T, root[1], 0, 1e-12)) { // 实根
        //                 const m_candidate = root[0];
        //                 if (2 * m_candidate - p >= -1e-12) { // 允许微小误差
        //                     m = m_candidate;
        //                     found_m = true;
        //                     break;
        //                 }
        //             }
        //         }
        //         if (!found_m) {
        //             // 无法找到合适的实 m，可能需要复数处理或返回失败
        //             @memset(re[0..4], [2]T{ 0, 0 });
        //             return re;
        //         }
        //         // 计算分解参数
        //         const s_squared = 2 * m - p;
        //         if (s_squared < 0) {
        //             @memset(@as([*][2]T, @ptrCast(re))[0..4], [2]T{ 0, 0 });
        //             return re;
        //         }

        //         const s = std.math.sqrt(s_squared);
        //         const t1 = (s * s + p - q / s) / 2;
        //         const t2 = (s * s + p + q / s) / 2;
        //         // 解两个二次方程
        //         var root_count: usize = 0;

        //         // 第一个二次方程: y^2 + s*y + t1 = 0
        //         const disc1 = s * s - 4 * t1;
        //         if (disc1 >= 0) {
        //             const sqrt_d = std.math.sqrt(disc1);
        //             re[root_count][0] = (-s + sqrt_d) / 2 - b_over_4;
        //             re[root_count][1] = 0;
        //             root_count += 1;
        //             if (root_count < 4 and disc1 > 0) {
        //                 re[root_count][0] = (-s - sqrt_d) / 2 - b_over_4;
        //                 re[root_count][1] = 0;
        //                 root_count += 1;
        //             }
        //         } else {
        //             const sqrt_d = std.math.sqrt(-disc1);
        //             re[root_count][0] = -s / 2 - b_over_4;
        //             re[root_count][1] = sqrt_d / 2;
        //             root_count += 1;
        //             re[root_count][0] = -s / 2 - b_over_4;
        //             re[root_count][1] = -sqrt_d / 2;
        //             root_count += 1;
        //         }

        //         // 第二个二次方程: y^2 - s*y + t2 = 0
        //         const disc2 = s * s - 4 * t2;
        //         if (disc2 >= 0) {
        //             const sqrt_disc2 = std.math.sqrt(disc2);
        //             // 两个实根
        //             if (root_count < 4) {
        //                 re[root_count][0] = (s + sqrt_disc2) / 2 - b_over_4;
        //                 re[root_count][1] = 0;
        //                 root_count += 1;
        //             }
        //             if (root_count < 4) {
        //                 re[root_count][0] = (s - sqrt_disc2) / 2 - b_over_4;
        //                 re[root_count][1] = 0;
        //                 root_count += 1;
        //             }
        //         } else {
        //             // 复根：y = [s ± i√|disc2|] / 2
        //             const sqrt_disc2 = std.math.sqrt(-disc2); // 正实数
        //             if (root_count < 4) {
        //                 re[root_count][0] = s / 2 - b_over_4; // 实部
        //                 re[root_count][1] = sqrt_disc2 / 2; // 虚部
        //                 root_count += 1;
        //             }
        //             if (root_count < 4) {
        //                 re[root_count][0] = s / 2 - b_over_4; // 实部
        //                 re[root_count][1] = -sqrt_disc2 / 2; // 虚部
        //                 root_count += 1;
        //             }
        //         }

        //         // 用0填充剩余位置
        //         while (root_count < 4) : (root_count += 1) {
        //             re[root_count][0] = 0;
        //             re[root_count][1] = 0;
        //         }
        //     }
        //     return re;
        // }
    };
}
